---
title: "Code and replication materials"
subtitle: "All in this together? A preregistered report on deservingness of government aid during the COVID-19 pandemic"
author: "Aengus Bridgman"
date: "March 2021"
output: pdf_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, cache = TRUE)
```

# Packages and functions

## Packages

```{r, echo = TRUE, results = 'hide', warning = FALSE, message = FALSE}

# renv::restore()
library(filesstrings)
library(cregg)
library(cjoint)
library(tidyverse)
library(haven)
library(patchwork)
library(texreg)
library(miceadds)
library(RColorBrewer)
library(scales)

# create outputs directory
dir.create("outputs/", showWarnings = FALSE)

# move data files if downloaded directly from dataverse
if (!file.exists("data/")) dir.create("data/")
if (file.exists("study1.csv")) file.rename("study1.csv", "data/study1.csv")
if (file.exists("study2.csv")) file.rename("study2.csv", "data/study2.csv")
if (file.exists("study2_alt.csv")) file.rename("study2_alt.csv", "data/study2_alt.csv")

`%nin%` <- Negate(`%in%`)

```

## Labelling functions

```{r, include = FALSE}

# These are recode functions for features and labels produced by the cregg:cj function and are used several times in the subsequent code. These are purely cosmetic.

recode_feature <- function(df) {
  
  df %>% 
    mutate(feature = case_when(
      str_detect(feature, "ethnicity") ~ "Ethnicity",
      str_detect(feature, "gender") ~ "Gender",
      str_detect(feature, "citizen") ~ "Citizenship\nstatus",
      str_detect(feature, "health") ~ "Health\nstatus",
      str_detect(feature, "marital") ~ "Marital\nstatus",
      str_detect(feature, "children") ~ "Children",
      str_detect(feature, "employ") ~ "Employment",
      str_detect(feature, "income") ~ "Previous\nincome"
    )) %>% 
    mutate(feature = factor(
      feature,
      levels = c("Previous\nincome","Health\nstatus","Employment",
                 "Children","Ethnicity","Citizenship\nstatus",
                 "Gender","Marital\nstatus")
    )) %>% 
    return()
}

relevel_study2 <- function(df) {
  df %>% 
    mutate(level = factor(
      level,
      levels = c("Female","Male",
                 "White","Visible minority",
                 "Healthy","Poor health", 
                 "Not a citizen","Citizen born in Canada","Citizen not born in Canada",
                 "Single","Married",
                 "No children","Children under 5",
                 "$120,000","$90,000","$60,000","$30,000",
                 "Employed full-time","Underemployed","Unemployed"
      ))) %>% 
    return()
}

refactor_amce <- function(df) {
  
  df %>% 
    mutate(feature = as.character(feature), level = as.character(level)) %>% 
    mutate(
      feature = case_when(
        str_detect(feature, "gender") ~ "Gender\n(ref.\nFemale)",
        str_detect(feature, "marital") ~ "Marital status\n(ref.\nSingle)",
        str_detect(feature, "health") ~ "Health\n(ref.\nHealthy)",
        str_detect(feature, "income") ~ "Income\n(ref.\n$120,000)",
        str_detect(feature, "ethnicity") ~ "Ethnicity\n(ref.\nWhite)",
        str_detect(feature, "employ") ~ "Employment\n(ref.\nEmployed FT)",
        str_detect(feature, "citizen") ~ "Citizenship\n(ref.\nNon-citizen)",
        str_detect(feature, "children") ~ "Children\n(ref.\nNo children)",
        TRUE ~ feature
    )) %>% 
      mutate(
        feature = factor(feature, 
                         levels = c("Income\n(ref.\n$120,000)",
                                    "Health\n(ref.\nHealthy)",
                                     "Employment\n(ref.\nEmployed FT)",
                                    "Children\n(ref.\nNo children)",
                                    "Ethnicity\n(ref.\nWhite)",
                                    "Citizenship\n(ref.\nNon-citizen)",
                                    "Gender\n(ref.\nFemale)",
                                    "Marital status\n(ref.\nSingle)"
                                    )),
        level = ordered(level,
                       levels = c("$90,000","$60,000","$30,000",
                                  "Citizen not born in Canada","Citizen born in Canada",
                                  "Children under 5",
                                  "Male","Female", "Single","Married", "Poor health",
                                  "Underemployed","Unemployed","Visible minority"))) %>% 
    return()
}


refactor_amce2 <- function(df) {
  
  df %>% 
    mutate(feature = as.character(feature), level = as.character(level)) %>% 
    mutate(
      feature = case_when(
        str_detect(feature, "gender") ~ "Gender\n(ref.\nFemale)",
        str_detect(feature, "marital") ~ "Marital status\n(ref.\nSingle)",
        str_detect(feature, "health") ~ "Health\n(ref.\nHealthy)",
        str_detect(feature, "income") ~ "Income\n(ref.\n$120,000)",
        str_detect(feature, "ethnicity") ~ "Ethnicity\n(ref.\nWhite)",
        str_detect(feature, "employ") ~ "Employment\n(ref.\nEmployed FT)",
        str_detect(feature, "citizen") ~ "Citizenship\n(ref.\nNon-citizen)",
        str_detect(feature, "children") ~ "Children\n(ref.\nNo children)",
        TRUE ~ feature,
    )) %>% 
      mutate(
        feature = factor(feature, 
                         levels = c("Income\n(ref.\n$120,000)",
                                    "Health\n(ref.\nHealthy)",
                                    "Employment\n(ref.\nEmployed FT)",
                                    "Children\n(ref.\nNo children)",
                                    "Ethnicity\n(ref.\nWhite)",
                                    "Citizenship\n(ref.\nNon-citizen)",
                                    "Gender\n(ref.\nFemale)",
                                    "Marital status\n(ref.\nSingle)"
                         ))) %>% 
    return()
}
```

## Theme

```{r}

palette_to_use <- brewer.pal(n = 11, name = "Spectral")

theme <- theme_bw() +
  theme(text = element_text(size = 12),
        axis.title = element_text(size = 12),
        plot.title = element_text(size = 12, face = "bold"),
        legend.position="top",
        legend.title = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.y = element_line(color = "lightgrey", size = 0.1, linetype = "solid"),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.x = element_line(color = "lightgrey", size = 0.2, linetype = "solid"),
        strip.background = element_blank(), strip.text = element_blank()
  ) 
```

## Data

# Study 1

## Figure 1: AMCEs and MMs

```{r, include = FALSE}

# Read in study1.csv, factor the id variable (which would otherwise be numeric) and several of the vignette_ variables. 
# Each observation is one id-vignette_num observation (id being an id number per respondent, and vignette_num being the order of the profiles shown).

study1 <- read.csv("data/study1.csv") %>% 
  mutate(
    id = factor(id),
    vignette_children = factor(
      vignette_children,
      levels = c("No children","Children under 5","Children over 12","Children between 5 and 12")
    ),
    vignette_citizen = factor(
      vignette_citizen,
      levels = c("Not a citizen","Citizen born in Canada","Citizen not born in Canada")
    ),
    vignette_employment = factor(
      vignette_employment,
      levels = c("Employed full-time","Employed, reduced income", "Unemployed due to pandemic","Unemployed")
    ),
    vignette_income = factor(
      vignette_income,
      levels = c("$120,000","$90,000","$60,000","$30,000")
    ),
    vignette_gender = factor(vignette_gender),
    vignette_health = factor(vignette_health),
    vignette_marital = factor(vignette_marital),
    vignette_ethnicity = factor(vignette_ethnicity)
  )

# Each Figure in the paper consists of 2-4 panels. The labelling is consistently figxx_y with xx refering to the figure number and y refering to the panel letter
fig1_A <-  study1 %>% 
  cregg::cj(id = ~id, estimate = "mm",
     formula = vignette_allocation ~ vignette_gender + vignette_citizen + vignette_health + vignette_marital + vignette_children + vignette_employment + vignette_income + vignette_ethnicity) %>% 
  recode_feature() %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  mutate(level = ordered(
    level,
    levels = c("$120,000","$90,000","$60,000","$30,000",
               "Healthy","Poor health",
               "Employed full-time","Employed, reduced income","Unemployed due to pandemic","Unemployed",
               "No children","Children under 5","Children between 5 and 12","Children over 12",
               "East Asian","White","South Asian","Indigenous",
               "Not a citizen","Citizen born in Canada","Citizen not born in Canada",
               "Female","Male",
               "Married","Single"))) %>% 
  ggplot(aes(x = level, y = estimate)) +
  geom_hline(aes(yintercept = mean(study1$vignette_allocation)), linetype = "dashed", color = "darkgrey") +
  geom_point() +
  scale_y_continuous(limits = c(1000,2000), breaks = seq(1000,2000, 200)) + 
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0) +
  theme +
  coord_flip() +
  scale_x_discrete(position = "top", labels = NULL) + 
  facet_grid(feature~., scales = "free", switch="both") +
  labs(x = "", y = "Allocation per month\nduring pandemic ($)",
       title = "(A) Marginal means") +
  theme(strip.background = element_blank(), axis.ticks = element_blank())
  
fig1_B <-  cjoint::amce(data = study1,
               formula = vignette_allocation ~ vignette_gender + vignette_citizen + vignette_health + vignette_marital + vignette_children + vignette_employment + vignette_income + vignette_ethnicity,
               na.ignore = TRUE, cluster = FALSE) %>% 
  .$estimates %>% 
  data.frame() %>% 
  t() %>% 
  data.frame() %>%
  rownames_to_column("attribute") %>% 
  # Split into feature and level using grep 
  mutate(feature = sub("?([A-Z]|[1-9]).*", "", attribute),
         level = sub("^.*?([A-Z]|[5])", "\\1", attribute)) %>% 
  select(-attribute) %>% 
  rownames_to_column() %>% 
  group_by(rowname) %>% 
  mutate(level = ifelse(is.na(level),
                        as.character(str_split(feature, "[.]", simplify = TRUE)[2]),
                        level),
         feature = as.character(str_split(feature, "[.]", n = 2, simplify = TRUE)[1])) %>% 
  ungroup() %>% 
  add_row(feature = "vignettechildren", level = "No children", AMCE = 0, Std..Error = 0) %>% 
  add_row(feature = "vignettecitizen", level = "Not a citizen", AMCE = 0, Std..Error = 0) %>% 
  add_row(feature = "vignetteemployment", level = "Employed full-time", AMCE = 0, Std..Error = 0) %>% 
  add_row(feature = "vignetteethnicity", level = "East Asian", AMCE = 0, Std..Error = 0) %>% 
  add_row(feature = "vignettegender", level = "Female", AMCE = 0, Std..Error = 0) %>% 
  add_row(feature = "vignettehealth", level = "Healthy", AMCE = 0, Std..Error = 0) %>% 
  add_row(feature = "vignetteincome", level = "$120,000", AMCE = 0, Std..Error = 0) %>% 
  add_row(feature = "vignettemarital", level = "Married", AMCE = 0, Std..Error = 0) %>% 
  recode_feature() %>% 
  mutate(level = ordered(case_when(
    str_detect(tolower(level), "5and12") ~ "Children between 5 and 12", 
    str_detect(tolower(level), "under5") ~ "Children under 5", 
    str_detect(tolower(level), "over12") ~ "Children over 12", 
    str_detect(tolower(level), "poorhealth") ~ "Poor health", 
    str_detect(tolower(level), "south") ~ "South Asian", 
    str_detect(tolower(level), "nochildren") ~ "No children", 
    str_detect(tolower(level), "citizennotbornin") ~ "Citizen not born in Canada", 
    str_detect(tolower(level), "citizenbornin") ~ "Citizen born in Canada", 
    str_detect(tolower(level), "unemployedduetopandemic") ~ "Unemployed due to pandemic", 
    str_detect(tolower(level), "reducedincome") ~ "Employed, reduced income", 
    str_detect(tolower(level), "90000") ~ "$90,000", 
    str_detect(tolower(level), "60000") ~ "$60,000", 
    str_detect(tolower(level), "30000") ~ "$30,000", 
    TRUE ~ level), 
    levels = c("$120,000","$90,000","$60,000","$30,000",
               "Healthy","Poor health",
               "Employed full-time","Employed, reduced income","Unemployed due to pandemic","Unemployed",
               "No children","Children under 5","Children between 5 and 12","Children over 12",
               "East Asian","White","South Asian","Indigenous",
               "Not a citizen","Citizen born in Canada","Citizen not born in Canada",
               "Female","Male",
               "Married","Single"))) %>% 
  mutate(lb = AMCE - 1.96*Std..Error,
         ub = AMCE + 1.96*Std..Error) %>% 
  ggplot(aes(x = level, y = AMCE)) +
  geom_hline(aes(yintercept = 0), linetype = "dashed", color = "darkgrey") +
  geom_point() +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0) +
  theme +
  scale_x_discrete(position = "top") + 
  coord_flip() +
  facet_grid(feature~., scales = "free", switch="both") +
  theme(strip.background = element_blank(), strip.text = element_blank()) +
  labs(title = "(B) AMCE", y = "AMCE of allocation per month\nduring pandemic ($)", x = "")

fig1_A + fig1_B

ggsave("outputs/fig1.pdf", width = 7.5, height = 7.5)
ggsave("outputs/fig1.png", width = 7.5, height = 7.5)

```

## Figure 2: Sub-group AMCEs 1

```{r, include = FALSE}

fig2_A <- bind_rows(
  filter(study1, respondent_employment == "At risk/laid off") %>% 
    cj(id = ~id, estimate = "mm",
       formula = vignette_allocation ~ vignette_employment,
       level_order = "descending") %>% 
    mutate(group = "At risk/laid off"),
  filter(study1, respondent_employment == "Not at risk") %>% 
    cj(id = ~id, estimate = "mm",
       formula = vignette_allocation ~ vignette_employment,
       level_order = "descending") %>% 
    mutate(group = "Not at risk")
) %>% 
  mutate(level = factor(level,
                        levels = c("Employed full-time","Employed, reduced income","Unemployed due to pandemic","Unemployed")),
         lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  ggplot(aes(x = level, y = estimate, col = group)) +
  geom_point(position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, position = position_dodge(0.3)) +
  coord_flip() +
  theme +
  theme(legend.title = element_blank(),
        legend.position="right") +
  labs(x = "", y = "") +
  guides(col = FALSE) +
  labs(title = "(A) Employment MMs") +
  scale_color_manual(values = palette_to_use[c(2,11)])

fig2_B <- bind_rows(
  study1 %>% 
    filter(respondent_employment == "At risk/laid off") %>% 
    cj(id = ~id, estimate = "amce",
       formula = vignette_allocation ~ vignette_employment,
       level_order = "descending") %>% 
    mutate(group = "At risk/laid off"),
  study1 %>% 
    filter(respondent_employment == "Not at risk") %>% 
    cj(id = ~id, estimate = "amce",
       formula = vignette_allocation ~ vignette_employment,
       level_order = "descending") %>% 
    mutate(group = "Not at risk")
) %>% 
  mutate(level = factor(level,
                        levels = c("Employed full-time","Employed, reduced income","Unemployed due to pandemic","Unemployed")),
         lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  ggplot(aes(x = level, y = estimate, col = group)) +
  geom_point(position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, position = position_dodge(0.3)) +
  theme +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.title = element_blank(),
        legend.position="right") +
  coord_flip() +
  labs(x = "", y = "") +
  labs(title = "(B) Employment AMCEs") +
  scale_color_manual(values = palette_to_use[c(2,11)])


fig2_C <- bind_rows(
  study1 %>% 
    filter(respondent_children == "No children") %>% 
    mutate(vignette_children = factor(ifelse(vignette_children == "No children", "No children", "Children"))) %>% 
    cj(id = ~id, estimate = "mm",
       formula = vignette_allocation ~ vignette_children,
       level_order = "descending") %>% 
    mutate(group = "No children"),
  study1 %>% 
    filter(respondent_children == "Children") %>% 
    mutate(vignette_children = factor(ifelse(vignette_children == "No children", "No children", "Children"))) %>% 
    cj(id = ~id, estimate = "mm",
       formula = vignette_allocation ~ vignette_children,
       level_order = "descending") %>% 
    mutate(group = "Children")
) %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  ggplot(aes(x = level, y = estimate, col = group)) +
  geom_point(position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, position = position_dodge(0.3)) +
  theme +
  coord_flip() +
  labs(x = "", y = "Allocation per month\nduring pandemic ($)") +
  guides(col = FALSE) +
  labs(title = "(C) Children MMs") +
  scale_color_manual(values = palette_to_use[c(2,11)])

fig2_D <- bind_rows(
  study1 %>% 
    filter(respondent_children == "Children") %>% 
    mutate(vignette_children = factor(ifelse(vignette_children == "No children", "No children", "Children"))) %>% 
    cj(id = ~id, estimate = "amce",
       formula = vignette_allocation ~ vignette_children,
       level_order = "descending") %>% 
    mutate(group = "No children"),
  study1 %>% 
    filter(respondent_children == "No children") %>% 
    mutate(vignette_children = factor(ifelse(vignette_children == "No children", "No children", "Children"))) %>% 
    cj(id = ~id, estimate = "amce",
       formula = vignette_allocation ~ vignette_children,
       level_order = "descending") %>% 
    mutate(group = "Children")
) %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  ggplot(aes(x = level, y = estimate, col = group)) +
  geom_point(position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, position = position_dodge(0.3)) +
  theme +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.title = element_blank()) +
  coord_flip() +
  labs(x = "", y = "AMCE of allocation per month\nduring pandemic ($)") +
  labs(title = "(D) Children AMCEs") +
  theme(legend.position="right") +
  scale_color_manual(values = palette_to_use[c(2,11)])

(fig2_A + fig2_B) / (fig2_C + fig2_D)

ggsave("outputs/fig2.pdf", width = 7.5, height = 5)
ggsave("outputs/fig2.png", width = 7.5, height = 5)

```

## Figure 3: Sub-group AMCEs 2

```{r, include = FALSE}


fig3_A <- bind_rows(
  filter(study1, respondent_comb_health == 0) %>% 
    cj(id = ~id, estimate = "mm", formula = vignette_allocation ~ vignette_health) %>% 
    mutate(group = "Good health"),
  filter(study1, respondent_comb_health == 1) %>% 
    cj(id = ~id, estimate = "mm", formula = vignette_allocation ~ vignette_health) %>% 
    mutate(group = "Poor Health")
) %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  ggplot(aes(x = level, y = estimate, col = group)) +
  geom_point(position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, position = position_dodge(0.3)) +
  theme +
  coord_flip() +
  guides(col = FALSE) +
  labs(x = "", y = "", title = "(A) Health MMs") +
  scale_color_manual(values = palette_to_use[c(2,11)])

fig3_B <- bind_rows(
    filter(study1, respondent_comb_health == 0) %>% 
    cj(id = ~id, estimate = "amce", formula = vignette_allocation ~ vignette_health) %>% 
    mutate(group = "Good health"),
    filter(study1, respondent_comb_health == 1) %>% 
    cj(id = ~id, estimate = "amce", formula = vignette_allocation ~ vignette_health) %>% 
    mutate(group = "Poor Health")
) %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  ggplot(aes(x = level, y = estimate, col = group)) +
  geom_point(position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, position = position_dodge(0.3)) +
  theme +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.title = element_blank(),
        legend.position="right") +
  coord_flip() +
  labs(x = "", y = "", title = "(B) Health AMCEs") +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.title = element_blank()) +
  scale_color_manual(values = palette_to_use[c(2,11)])

fig3_C <- bind_rows(
    filter(study1, respondent_income == "Low income") %>% 
    cj(id = ~id, estimate = "mm", 
       formula = vignette_allocation ~ vignette_income,
       level_order = "descending") %>% 
    mutate(group = "Low income"),
    filter(study1, respondent_income == "High income") %>% 
    cj(id = ~id, estimate = "mm", 
       formula = vignette_allocation ~ vignette_income,
       level_order = "descending") %>% 
    mutate(group = "High income")
) %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  ggplot(aes(x = level, y = estimate, col = group)) +
  geom_point(position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, position = position_dodge(0.3)) +
  theme +
  guides(col = FALSE) +
  coord_flip() +
  labs(x = "", y = "Allocation per month\nduring pandemic ($)", title = "(C) Income MMs") +
  scale_color_manual(values = palette_to_use[c(2,11)])

fig3_D <- bind_rows(
  filter(study1, respondent_income == "Low income") %>% 
    cj(id = ~id, estimate = "amce", 
       formula = vignette_allocation ~ vignette_income,
       level_order = "descending") %>% 
    mutate(group = "Low income"),
  filter(study1, respondent_income == "High income") %>% 
    cj(id = ~id, estimate = "amce",        
       formula = vignette_allocation ~ vignette_income,
       level_order = "descending") %>% 
    mutate(group = "High income")
) %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  ggplot(aes(x = level, y = estimate, col = group)) +
  geom_point(position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, position = position_dodge(0.3)) +
  theme +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.title = element_blank(),
        legend.position="right") +
  coord_flip() +
  labs(x = "", y = "AMCE of allocation per month\nduring pandemic ($)", title = "(D) Income AMCEs") +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.title = element_blank()) +
  scale_color_manual(values = palette_to_use[c(2,11)])


(fig3_A + fig3_B) / (fig3_C + fig3_D)

ggsave("outputs/fig3.pdf", width = 7.5, height = 5)
ggsave("outputs/fig3.png", width = 7.5, height = 5)

```

# Study 2

```{r}

study2 <- read.csv("data/study2.csv") %>% 
  group_by(id) %>% 
  mutate(
    vignette_gender = factor(vignette_gender, levels = c("Female","Male")),
    vignette_ethnicity = factor(vignette_ethnicity, levels = c("White","Visible minority")),
    vignette_health = factor(vignette_health, levels = c("Healthy","Poor health")),
    vignette_citizen = factor(vignette_citizen, levels = c("Not a citizen","Citizen born in Canada",
                                                           "Citizen not born in Canada")),
    vignette_marital = factor(vignette_marital, levels = c("Single","Married")),
    vignette_children = factor(vignette_children, levels = c("No children","Children under 5")),
    vignette_income = factor(vignette_income, levels = c("$120,000","$90,000","$60,000","$30,000")),
    vignette_employment = factor(vignette_employment, levels = c("Employed full-time","Underemployed","Unemployed")),
    
    # Create two new variables meaning over the responses around support for government redistribution 
    gov_covid = mean(c(gov_cerb,gov_cews,gov_cesb,gov_ccb), na.rm = TRUE),
    gov_overall = mean(c(gov_cerb,gov_cews,gov_cesb,gov_tax,gov_ccb,gov_ei,gov_support), na.rm = TRUE)) %>% 
  ungroup() %>% 
  mutate(id = factor(id))


# Remove respondents below a third of the median time and who replied $9999 or above to a monthly allocation amount
median_completion <- median(study2$duration)
study2 <- filter(study2, duration > median_completion*0.33)
study2 <- filter(study2, vignette_allocation < 9999)

# Also exclude the top 2 percent
exclusion_criteria <- study2 %>% 
  group_by(condition) %>%
  summarize(exclude = quantile(vignette_allocation, 0.98, na.rm = TRUE), .groups = "drop")
max_covid_allocation <- as.numeric(exclusion_criteria[1,2])
max_gst_allocation <- as.numeric(exclusion_criteria[2,2])

study2 <- filter(study2,
                 condition == "COVID-19" & vignette_allocation < max_covid_allocation |
                   condition == "GST" & vignette_allocation < max_gst_allocation)

# Show descriptive statistics for each condition
study2 %>% 
  group_by(condition) %>%
  summarize(n = n(), 
            sd = sd(vignette_allocation, na.rm = TRUE),
            mean = mean(vignette_allocation, na.rm = TRUE),
            median = median(vignette_allocation, na.rm = TRUE),
            max = max(vignette_allocation, na.rm = TRUE),
            .groups = "drop")

```

## Figure 4: AMCEs

```{r}

# Get the mean for each condition to use to normalize the data
mean_covid <- mean(filter(study2, condition == "COVID-19")$vignette_allocation)
mean_gst<- mean(filter(study2, condition == "GST")$vignette_allocation)

study2 <- study2 %>% 
  mutate(allocation_normalized = case_when(
    condition == "COVID-19" ~ (vignette_allocation/mean_covid),
    condition == "GST" ~ (vignette_allocation/mean_gst)
  )) 

fig4_A <- study2 %>% 
  filter(condition == "COVID-19") %>% 
  cj(., id = ~id, estimate = "amce",
     vignette_allocation ~ vignette_gender + vignette_ethnicity + vignette_health + vignette_citizen + vignette_marital + vignette_children + vignette_employment + vignette_income) %>% 
  recode_feature() %>% 
  relevel_study2() %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  ggplot(aes(x = level, y = estimate), color = "black", fill = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") + 
  geom_point(color = palette_to_use[3]) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, color = palette_to_use[3]) +
  theme +
  coord_flip() +
  scale_x_discrete(position = "top", labels = NULL) + 
  scale_y_continuous(limits = c(-0.2*mean_covid, mean_covid), breaks = c(-200, 0, 400, 800, 1200)) +
  facet_grid(feature~., scales = "free", switch="both") +
  labs(x = "", y = "AMCE ($), with axis\nlimit set at mean of\nCOVID-19 allocation",
       title = "(A) COVID-19 allocation") +
  theme + 
  theme(strip.background = element_blank(), axis.ticks = element_blank())

fig4_B <- study2 %>% 
  filter(condition == "GST") %>% 
  cj(., id = ~id, estimate = "amce",
     vignette_allocation ~ vignette_gender + vignette_ethnicity + vignette_health + vignette_citizen + vignette_marital + vignette_children + vignette_employment + vignette_income) %>% 
  recode_feature() %>% 
  relevel_study2() %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  ggplot(aes(x = level, y = estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") + 
  geom_point(color = palette_to_use[10]) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, color = palette_to_use[10]) +
  theme +
  coord_flip() +
  scale_x_discrete(position = "top") + 
  #scale_y_continuous(limits = c(-0.2*me_gst, me_gst)) +
  facet_grid(feature~., scales = "free", switch="both") +
  labs(x = "", y = "AMCE ($), with axis\nlimit set at mean of\nGST allocation", title = "(B) GST rebate") +
  theme +
  theme(strip.text = element_blank(), strip.background = element_blank(), axis.ticks = element_blank())

fig4_A + fig4_B

ggsave("outputs/fig4.pdf", width = 7.5, height = 7.5)
ggsave("outputs/fig4.png", width = 7.5, height = 7.5)

```

## Figure 5: AMCE differences

```{r}

fig5_A <- bind_rows(
  filter(study2, condition == "COVID-19") %>% 
    cj(., id = ~id, estimate = "amce", level_order = "descending",
         allocation_normalized ~ vignette_gender + vignette_ethnicity + vignette_health + vignette_citizen + vignette_marital + vignette_children + vignette_employment + vignette_income) %>% 
    mutate(condition = "COVID-19"),
  filter(study2, condition == "GST") %>% 
    cj(., id = ~id, estimate = "amce", level_order = "descending",
         allocation_normalized ~ vignette_gender + vignette_ethnicity + vignette_health + vignette_citizen + vignette_marital + vignette_children + vignette_employment + vignette_income) %>% 
    mutate(condition = "GST")
) %>% 
  mutate(feature = as.character(feature)) %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  filter(estimate != 0) %>% 
  refactor_amce() %>% 
  ggplot(aes(x = level, y = estimate, fill = condition, col = condition)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
  geom_point(position = position_dodge(0.8)) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, position = position_dodge(0.8)) +
  theme +
  coord_flip() +
  scale_color_manual(values = palette_to_use[c(3,10)]) + 
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_x_discrete(position = "top", labels = NULL) + 
  facet_grid(feature~., scales = "free", switch="both") +
  labs(x = "", y = "AMCE (% mean)", title = "(A) AMCEs") +
  theme(strip.background = element_blank(), 
        axis.ticks = element_blank(),
        strip.text = element_text(size = 8),
        legend.title = element_blank(), legend.position = "bottom")


fig5_B <- study2 %>% 
  amce_diffs(., id = ~id, by = ~condition, level_order = "descending",
     allocation_normalized ~ vignette_gender + vignette_ethnicity + vignette_health + vignette_citizen + vignette_marital + vignette_children + vignette_employment + vignette_income) %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  refactor_amce() %>% 
  ggplot(aes(x = level, y = estimate, col = "Difference")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
  geom_point(col = "black") +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, col = "black") +
  theme +
  coord_flip() +
  scale_x_discrete(position = "top") + 
  scale_color_manual(values = "black") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + 
  facet_grid(feature~., scales = "free", switch="both") +
  labs(x = "", y = "AMCE difference (% mean,\nGST minus COVID-19)", title = "(B) AMCE differences") +
  theme(strip.background = element_blank(), 
        strip.text = element_blank(), axis.ticks = element_blank()) + theme(legend.position = "bottom", legend.title = element_blank())

fig5_A + fig5_B

ggsave("outputs/fig5.pdf", width = 7.5, height = 8.5)
ggsave("outputs/fig5.png", width = 7.5, height = 8.5)

```

## Figure 6: Deservingness and similarity

```{r}

mean_covid <- mean(filter(study2, condition == "COVID-19")$vignette_allocation)
mean_gst<- mean(filter(study2, condition == "GST")$vignette_allocation)

study2 <- mutate(study2, 
                 allocation_normalized = vignette_allocation/ifelse(condition == "COVID-19", mean_covid, mean_gst))


set.seed(0)
fig6_A <- filter(study2, condition == "COVID-19") %>% 
  ggplot(aes(x = vignette_deservingness, y = allocation_normalized, fill = condition)) +
  geom_jitter(col = palette_to_use[c(3)], alpha = 0.3, width = 0.25, height = 0.25) +
  geom_smooth(formula = y ~ x, method='lm', color = "black", size = 0.75, fill = palette_to_use[c(3)], alpha = 1) +
  scale_x_continuous(breaks = seq(0,10,1)) +
  labs(x = "", y = "Allocation\n(% of mean, $1152)", title = "(A) COVID-19 deservingness") +
  theme + 
  theme(strip.background = element_blank()) + 
  scale_y_continuous(labels = scales::percent) +
  coord_cartesian(ylim = c(0, 3.5)) + 
  guides(fill = FALSE, color = FALSE)

set.seed(0)
fig6_B <- filter(study2, condition == "COVID-19") %>% 
  ggplot(aes(x = vignette_similarity, y = allocation_normalized, fill = condition)) +
  geom_jitter(col = palette_to_use[3], alpha = 0.3, width = 0.25, height = 0.25) +
  geom_smooth(formula = y ~ x, method='lm', color = "black", size = 0.75, fill = palette_to_use[c(3)], alpha = 1) +
  scale_x_continuous(breaks = seq(0,10,1)) +
  labs(x = "", y = "", title = "(B) COVID-19 similarity") +
  theme + 
  theme(strip.background = element_blank()) + 
  scale_y_continuous(labels = scales::percent) +
  coord_cartesian(ylim = c(0, 3.5)) + 
  guides(fill = FALSE, color = FALSE)

set.seed(0)
fig6_C <- filter(study2, condition == "GST") %>% 
  ggplot(aes(x = vignette_deservingness, y = allocation_normalized, fill = condition)) +
  geom_jitter(col = palette_to_use[c(10)], alpha = 0.3, width = 0.25, height = 0.25) +
  geom_smooth(formula = y ~ x, method='lm', color = "black", size = 0.75, fill = palette_to_use[c(10)], alpha = 1) +
  scale_x_continuous(breaks = seq(0,10,1)) +
  labs(x = "Subjective deservingness", y = "Allocation\n(% of mean, $198)", title = "(C) GST deservingness") +
  theme + 
  theme(strip.background = element_blank()) + 
  scale_y_continuous(labels = scales::percent) +
  coord_cartesian(ylim = c(0, 3.5)) + 
  guides(fill = FALSE, color = FALSE)

set.seed(0)
fig6_D <- filter(study2, condition == "GST") %>% 
  ggplot(aes(x = vignette_similarity, y = allocation_normalized, fill = condition)) +
  geom_jitter(col = palette_to_use[10], alpha = 0.3, width = 0.25, height = 0.25) +
  geom_smooth(formula = y ~ x, method='lm', color = "black", size = 0.75, fill = palette_to_use[c(10)], alpha = 1) +
  scale_x_continuous(breaks = seq(0,10,1)) +
  labs(x = "Subjective similarity", y = "", title = "(D) GST similarity") +
  theme + 
  theme(strip.background = element_blank()) + 
  scale_y_continuous(labels = scales::percent) +
  coord_cartesian(ylim = c(0, 3.5)) + 
  guides(fill = FALSE, color = FALSE)


(fig6_A + fig6_B) / (fig6_C + fig6_D)

ggsave("outputs/fig6.pdf", width = 7.5, height = 4)
ggsave("outputs/fig6.png", width = 7.5, height = 4)

```

## F-test

```{r}

cj_anova(study2,
         id = ~id, by = ~condition,
         formula = allocation_normalized ~ vignette_ethnicity + vignette_health + 
           vignette_citizen + vignette_children + vignette_employment + vignette_income)


```

## Table 2: Subjective

```{r}

tab2_models <- list()

# Clustered standard errors at the individual level
tab2_models[[1]] <-  miceadds::lm.cluster(data = filter(study2, condition == "COVID-19"),
                  vignette_allocation ~ vignette_deservingness + vignette_similarity + as.factor(id), cluster = "id")
tab2_models[[2]] <- miceadds::lm.cluster(data = filter(study2, condition == "GST"),
                  vignette_allocation ~ vignette_deservingness + vignette_similarity + as.factor(id), cluster = "id") 
tab2_models[[3]] <- miceadds::lm.cluster(data = study2,
                  allocation_normalized ~ condition*vignette_deservingness + condition*vignette_similarity + as.factor(id), cluster = "id") 

texreg(list(tab2_models[[1]]$lm_res, 
                    tab2_models[[2]]$lm_res,
                    tab2_models[[3]]$lm_res),
               omit.coef = "id",
               file = "outputs/tab2.tex",
               single.row = TRUE,
               custom.model.names = c("COVID-19",
                                      "GST rebate",
                                      "Mean-normalized"),
               custom.coef.names = c("Constant",
                                     "Deservingness (0-10)",
                                     "Similarity (0-10)",
                                     "GST",
                                     "GST x Deservingness (0-10)",
                                     "GST x Similarity (0-10)"),
               reorder.coef = c(1,4,2,3,5,6),
               dcolumn = FALSE,
               stars = c(0.01),
               label = "tab:2",
               caption = "Subjective evaluations of deservingness and similarity",
               caption.above = TRUE,
               custom.note = "\\parbox{0.9\\linewidth}{\\vspace{2pt}%stars. Linear regression for subjective evaluations of deservingness and similarity, with individual respondent controls and clustered standard errors in parentheses. Dependent variable: allocation of cash transfer to hypothetical individuals, either under the COVID-19 or GST conditions.}")
```

# Appendix

## Appendix A: Study 1 objective similarity measures

### Figure A1: Sub-group AMCEs 3

```{r}

figA1_A <- bind_rows(
    filter(study1, respondent_gender == "Female") %>% 
    cj(id = ~id, estimate = "mm",
       formula = vignette_allocation ~ vignette_gender,
       level_order = "descending") %>% 
    mutate(group = "Female"),
  filter(study1, respondent_gender == "Male") %>% 
    cj(id = ~id, estimate = "mm",
       formula = vignette_allocation ~ vignette_gender,
       level_order = "descending") %>% 
    mutate(group = "Male")
) %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  ggplot(aes(x = level, y = estimate, col = group)) +
  geom_point(position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, position = position_dodge(0.3)) +
  theme +
  coord_flip() +
  labs(x = "", y = "", title = "(A) Gender MMs") +
  guides(col = FALSE) +
  scale_color_manual(values = palette_to_use[c(2,11)])

figA1_B <- bind_rows(
  filter(study1, respondent_gender == "Female") %>% 
    cj(id = ~id, estimate = "amce",
       formula = vignette_allocation ~ vignette_gender,
       level_order = "descending") %>% 
    mutate(group = "Female"),
  filter(study1, respondent_gender == "Male") %>% 
    cj(id = ~id, estimate = "amce",
       formula = vignette_allocation ~ vignette_gender,
       level_order = "descending") %>% 
    mutate(group = "Male")
) %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  ggplot(aes(x = level, y = estimate, col = group)) +
  geom_point(position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, position = position_dodge(0.3)) +
  theme +
  theme(legend.position="right") +
  coord_flip() +
  labs(x = "", y = "", title = "(B) Gender AMCEs") +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.title = element_blank()) +
  scale_color_manual(values = palette_to_use[c(2,11)])

figA1_C <- bind_rows(
  filter(study1, respondent_marital == "Married") %>% 
    cj(id = ~id, estimate = "mm",
       formula = vignette_allocation ~ vignette_marital,
       level_order = "descending") %>% 
    mutate(group = "Married"),
  filter(study1, respondent_marital == "Not married") %>% 
    cj(id = ~id, estimate = "mm",
       formula = vignette_allocation ~ vignette_marital,
       level_order = "descending") %>% 
    mutate(group = "Not married")
) %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  ggplot(aes(x = level, y = estimate, col = group)) +
  geom_point(position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, position = position_dodge(0.3)) +
  theme +
  coord_flip() +
  guides(col = FALSE) +
  labs(x = "", y = "Allocation per month\nduring pandemic ($)", title = "(C) Marital MMs") +
  scale_color_manual(values = palette_to_use[c(2,11)])

figA1_D <- bind_rows(
  filter(study1, respondent_marital == "Married") %>% 
    cj(id = ~id, estimate = "amce",
       formula = vignette_allocation ~ vignette_marital,
       level_order = "descending") %>% 
    mutate(group = "Married"),
  filter(study1, respondent_marital == "Not married") %>% 
    cj(id = ~id, estimate = "amce",
       formula = vignette_allocation ~ vignette_marital,
       level_order = "descending") %>% 
    mutate(group = "Not married")
) %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  ggplot(aes(x = level, y = estimate, col = group)) +
  geom_point(position = position_dodge(0.3)) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, position = position_dodge(0.3)) +
  theme +
  theme(legend.position="right") +
  coord_flip() +
  labs(x = "", y = "Allocation AMCEs per month\nduring pandemic ($)", title = "(D) Marital AMCEs") +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.title = element_blank()) +
  scale_color_manual(values = palette_to_use[c(2,11)])

(figA1_A + figA1_B) / (figA1_C + figA1_D)

ggsave("outputs/figA1.pdf", width = 7.5, height = 5)
ggsave("outputs/figA1.png", width = 7.5, height = 5)

```

## Appendix B: General support for redistribution

### Figure A2: outcomes and overall level of support for redistribution

```{r}

set.seed(0)
figA2_A <- study2 %>%
  filter(!is.na(gov_covid)) %>% 
  ggplot(aes(x = gov_covid, y = vignette_deservingness)) +
  geom_jitter(alpha = 0.05, width = 0.5, height = 0.5) + 
  geom_smooth(formula = y ~ x, method='lm', size = 0.75, color = "black", fill = "lightgrey") +
  theme +
  scale_y_continuous(breaks = seq(0,10,1)) +
  labs(x = "", y = "Subjective deservingness", title = "(A) Deservingness")

set.seed(0)
figA2_B <- study2 %>%
  filter(!is.na(gov_covid)) %>% 
  ggplot(aes(x = gov_covid, y = vignette_similarity)) +
  geom_jitter(alpha = 0.05, width = 0.5, height = 0.5) +
  geom_smooth(formula = y ~ x, method='lm', size = 0.75, color = "black", fill = "lightgrey") +
  theme +
  scale_y_continuous(breaks = seq(0,10,1)) +
  labs(x = "", y = "Subjective similarity", title = "(B) Similarity")

set.seed(0)
figA2_C <- study2 %>% 
  filter(condition == "COVID-19", !is.na(gov_covid)) %>%
  ggplot(aes(x = gov_covid, y = allocation_normalized, col = condition)) +
  geom_jitter(alpha = 0.1, width = 0.5, height = 0.3) +
  geom_smooth(formula = y ~ x, method='lm', size = 0.75, color = "black", fill = "lightgrey") +
  theme +
  scale_color_manual(values = palette_to_use[3]) + 
  scale_y_continuous(labels = scales::percent) +
  coord_cartesian(ylim = c(0,3)) + 
  labs(x = "Support for COVID-19 related transfers", y = "Normalized allocation amount\n(as percent of mean)", title = "(C) COVID-19") +
  guides(color =FALSE)

set.seed(0)
figA2_D <- study2 %>% 
  filter(condition == "GST", !is.na(gov_covid)) %>%
  ggplot(aes(x = gov_covid, y = allocation_normalized, col = condition)) +
  geom_jitter(alpha = 0.1, width = 0.5, height = 0.3) +
  geom_smooth(formula = y ~ x, method='lm', size = 0.75, color = "black", fill = "lightgrey") +
  theme +
  scale_color_manual(values = palette_to_use[10]) + 
  scale_y_continuous(labels = scales::percent) +
  coord_cartesian(ylim = c(0,3)) + 
  labs(x = "Support for COVID-19 related transfers", y = "Normalized allocation amount\n(as percent of mean)", title = "(D) GST") +
  guides(color =FALSE)


(figA2_A + figA2_B) / (figA2_C + figA2_D)

ggsave("outputs/figA2.pdf", width = 7.5, height = 7)
ggsave("outputs/figA2.png", width = 7.5, height = 7)

```

### Table A2: allocation and general support

```{r}

tabA2_models <- list()

tabA2_models[[1]] <- lm(data = study2,
                       allocation_normalized ~ gov_cerb)
tabA2_models[[2]] <- lm(data = study2,
                       allocation_normalized ~ gov_cews)
tabA2_models[[3]] <- lm(data = study2,
                       allocation_normalized ~ gov_support)
tabA2_models[[4]] <- lm(data = study2,
                       allocation_normalized ~ gov_ei)

texreg(tabA2_models,
               file = "outputs/tabA2.tex",
               single.row = TRUE,
               custom.coef.names = c("Constant",
                                     "Canadian Emergency Response Benefit",
                                     "Canadian Emergency Wage Subsidy",
                                     "General support for universal cash transfers",
                                     "General support for Employment Insurance enhancements"),
               dcolumn = FALSE,
               stars = c(0.01),
               label = "tab:a2",
               caption = "Correlations between mean-normalized allocations and general support for government redistribution programs",
               caption.above = TRUE,
               custom.note = "\\parbox{1\\linewidth}{\\vspace{2pt}%stars. Dependent variable: mean-normalized allocation of cash transfer to hypothetical individuals, either under the COVID-19 or GST conditions}")
```

### Table A3: allocation and attitudes about redistribution

```{r}

tabA3_models <- list()

tabA3_models[[1]] <- lm(data = study2,
                       allocation_normalized ~ get_ahead)
tabA3_models[[2]] <- lm(data = study2,
                       allocation_normalized ~ inequality)
tabA3_models[[3]] <- lm(data = study2,
                       allocation_normalized ~ living)

texreg(tabA3_models,
               file = "outputs/tabA3.tex",
               single.row = TRUE,
               custom.coef.names = c("Constant",
                                     "People who don't get ahead should blame themselves, not the system",
                                     "Government should take measures to reduce differences in income levels",
                                     "The government should see to it that everyone has a decent standard of living"),
               dcolumn = FALSE,
               stars = c(0.01),
               label = "tab:a3",
               caption = "Correlations between mean-normalized allocations and attitudes about redistribution",
               caption.above = TRUE,
               custom.note = "\\parbox{0.8\\linewidth}{\\vspace{2pt}%stars. Dependent variable: mean-normalized allocation of cash transfer to hypothetical individuals, either under the COVID-19 or GST conditions.}")


```

## Appendix C: Pre-peer reviewed experiment

```{r}

# This code is an exact replicate of that for Study 2

study2_alt <- read.csv("data/study2_alt.csv") %>% 
      mutate(
        vignette_gender = factor(vignette_gender, levels = c("Female","Male")),
        vignette_ethnicity = factor(vignette_ethnicity, levels = c("White","Visible minority")),
        vignette_health = factor(vignette_health, levels = c("Healthy","Poor health")),
        vignette_citizen = factor(vignette_citizen, levels = c("Not a citizen","Citizen born in Canada",
                                                               "Citizen not born in Canada")),
        vignette_marital = factor(vignette_marital, levels = c("Single","Married")),
        vignette_children = factor(vignette_children, levels = c("No children","Children under 5")),
        vignette_income = factor(vignette_income, levels = c("$120,000","$90,000","$60,000","$30,000")),
        vignette_employment = factor(vignette_employment, levels = c("Employed full-time","Underemployed","Unemployed"))
      ) %>% 
  mutate(id = factor(id))

median_completion <- median(study2_alt$duration)

study2_alt <- filter(study2_alt, duration > median_completion*0.33); nrow(study2_alt)
study2_alt <- filter(study2_alt, vignette_allocation < 9999); nrow(study2_alt)

# Exclusion
exclusion_criteria <- study2_alt %>% 
  group_by(condition) %>%
  summarize(exclude = quantile(vignette_allocation, 0.98, na.rm = TRUE), .groups = "drop")

max_covid_allocation <- as.numeric(exclusion_criteria[1,2])
max_gst_allocation <- as.numeric(exclusion_criteria[2,2])

study2_alt <- filter(study2_alt,
                 condition == "COVID-19" & vignette_allocation < max_covid_allocation |
                   condition == "GST" & vignette_allocation < max_gst_allocation)

study2_alt %>% 
  group_by(condition) %>%
  summarize(n = n(), 
            sd = sd(vignette_allocation, na.rm = TRUE),
            mean = mean(vignette_allocation, na.rm = TRUE),
            median = median(vignette_allocation, na.rm = TRUE),
            max = max(vignette_allocation, na.rm = TRUE),
            .groups = "drop")

median_completion <- median(study2_alt$duration)

study2_alt <- filter(study2_alt, duration > median_completion*0.33)
study2_alt <- filter(study2_alt, vignette_allocation < 9999); nrow(study2_alt)

# Exclusion 1
exclusion_criteria <- study2_alt %>% 
  group_by(condition) %>%
  summarize(exclude = quantile(vignette_allocation, 0.98, na.rm = TRUE), .groups = "drop")

max_covid_allocation <- as.numeric(exclusion_criteria[1,2])
max_gst_allocation <- as.numeric(exclusion_criteria[2,2])

study2_alt <- filter(study2_alt,
             condition == "COVID-19" & vignette_allocation < max_covid_allocation |
               condition == "GST" & vignette_allocation < max_gst_allocation)

study2_alt %>% 
  group_by(condition) %>%
  summarize(n = n(), 
            sd = sd(vignette_allocation, na.rm = TRUE),
            mean = mean(vignette_allocation, na.rm = TRUE),
            median = median(vignette_allocation, na.rm = TRUE),
            max = max(vignette_allocation, na.rm = TRUE),
            .groups = "drop")

```

### Figure A3: AMCEs

```{r}

figA3_A <- study2_alt %>% 
  filter(condition == "COVID-19") %>% 
  cj(., id = ~id, estimate = "amce",
     vignette_allocation ~ vignette_gender + vignette_ethnicity + vignette_health + vignette_citizen + vignette_marital + vignette_children + vignette_employment + vignette_income) %>% 
  recode_feature() %>% 
  relevel_study2() %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  ggplot(aes(x = level, y = estimate), color = "black", fill = "black") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
  geom_point(color = "black", fill = "black") +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, color = "black") +
  theme +
  coord_flip() +
  scale_x_discrete(position = "top", labels = NULL) + 
  facet_grid(feature~., scales = "free", switch="both") +
  labs(x = "", y = "AMCE of allocation per month\nduring pandemic ($)", title = "(A) COVID-19 allocation") +
  theme + 
  theme(strip.background = element_blank(), axis.ticks = element_blank())

figA3_B <- study2_alt %>% 
  filter(condition == "GST") %>% 
  cj(., id = ~id, estimate = "amce",
     vignette_allocation ~ vignette_gender + vignette_ethnicity + vignette_health + vignette_citizen + vignette_marital + vignette_children + vignette_employment + vignette_income) %>% 
  recode_feature() %>% 
  relevel_study2() %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  ggplot(aes(x = level, y = estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") + 
  geom_point(color = "black", fill = "black") +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, color = "black") +
  coord_flip() +
  scale_x_discrete(position = "top") + 
  facet_grid(feature~., scales = "free", switch="both") +
  labs(x = "", y = "AMCE of allocation per month\nduring pandemic ($)", title = "(B) GST rebate") +
  theme +
  theme(strip.text = element_blank(), strip.background = element_blank(), axis.ticks = element_blank())

figA3_A + figA3_B

ggsave("outputs/figA3.pdf", width = 7.5, height = 7.5)
ggsave("outputs/figA3.png", width = 7.5, height = 7.5)

```

### FigureA4: AMCE differences

```{r}

me_cov <- filter(study2_alt, condition == "COVID-19") %>% pull(vignette_allocation) %>% mean()
me_gst <- filter(study2_alt, condition == "GST") %>% pull(vignette_allocation) %>% mean()

study2_alt <- study2_alt %>% 
  mutate(allocation_normalized = case_when(
    condition == "COVID-19" ~ (vignette_allocation/me_cov),
    condition == "GST" ~ (vignette_allocation/me_gst)
  )) 

figA4_A <- bind_rows(
  filter(study2_alt, condition == "COVID-19") %>% 
    cj(., id = ~id, estimate = "amce", level_order = "descending",
         allocation_normalized ~ vignette_gender + vignette_ethnicity + vignette_health + vignette_citizen + vignette_marital + vignette_children + vignette_employment + vignette_income) %>% 
    mutate(condition = "COVID-19"),
  filter(study2_alt, condition == "GST") %>% 
    cj(., id = ~id, estimate = "amce", level_order = "descending",
         allocation_normalized ~ vignette_gender + vignette_ethnicity + vignette_health + vignette_citizen + vignette_marital + vignette_children + vignette_employment + vignette_income) %>% 
    mutate(condition = "GST")
) %>% 
  mutate(feature = as.character(feature)) %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  filter(estimate != 0) %>% 
  refactor_amce2() %>% 
  relevel_study2() %>% 
  ggplot(aes(x = level, y = estimate, fill = condition, col = condition)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
  geom_point(position = position_dodge(0.8)) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, position = position_dodge(0.8)) +
  theme +
  coord_flip() +
  scale_color_manual(values = palette_to_use[c(3,10)]) + 
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_x_discrete(position = "top", labels = NULL) + 
  facet_grid(feature~., scales = "free", switch="both") +
  labs(x = "", y = "AMCE (% mean allocation)", title = "(A) AMCE") +
  theme(strip.background = element_blank(), 
        axis.ticks = element_blank(),
        strip.text = element_text(size = 8),
        legend.title = element_blank(),
        legend.position = "bottom")

figA4_B <- study2_alt %>% 
  amce_diffs(., id = ~id, by = ~condition, level_order = "descending",
         allocation_normalized ~ vignette_gender + vignette_ethnicity + vignette_health + vignette_citizen + vignette_marital + vignette_children + vignette_employment + vignette_income) %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  refactor_amce2() %>% 
  relevel_study2() %>% 
  ggplot(aes(x = level, y = estimate, col = "Difference")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
  geom_point(col = "black") +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, col = "black") +
  theme +
  coord_flip() +
  scale_x_discrete(position = "top") + 
  scale_color_manual(values = "black") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + 
  facet_grid(feature~., scales = "free", switch="both") +
  labs(x = "", y = "AMCE difference (% mean allocation)", title = "(B) AMCE difference") +
  theme(strip.background = element_blank(), strip.text = element_blank(), axis.ticks = element_blank()) + theme(legend.position = "bottom", legend.title = element_blank())

figA4_A + figA4_B

ggsave("outputs/figA4.pdf", width = 7.5, height = 8.5)
ggsave("outputs/figA4.png", width = 7.5, height = 8.5)

```

### Figure A5: Deservingness and similarity

```{r}


mean_covid <- mean(filter(study2_alt, condition == "COVID-19")$vignette_allocation)
mean_gst<- mean(filter(study2_alt, condition == "GST")$vignette_allocation)

study2_alt <- mutate(study2_alt, allocation_normalized = vignette_allocation/ifelse(condition == "COVID-19", mean_covid, mean_gst))

mean_covid <- mean(filter(study2_alt, condition == "COVID-19")$vignette_allocation)
mean_gst<- mean(filter(study2_alt, condition == "GST")$vignette_allocation)

set.seed(0)
figA5_A <- filter(study2_alt, condition == "COVID-19") %>% 
  ggplot(aes(x = vignette_deservingness, y = allocation_normalized, fill = condition) ) +
  geom_jitter(col = palette_to_use[c(3)], alpha = 0.3, width = 0.25, height = 0.25) +
  geom_smooth(formula = y ~ x, method='lm', color = "black", size = 0.75, fill = palette_to_use[3], alpha = 1) +
  scale_x_continuous(breaks = seq(0,10,1)) +
  labs(x = "", y = "Allocation\n(% of mean, $1152)", title = "(A) COVID-19 deservingness") +
  theme + 
  theme(strip.background = element_blank()) + 
  scale_y_continuous(labels = scales::percent) +
  coord_cartesian(ylim = c(0, 3.5)) + 
  guides(fill = FALSE, color = FALSE)

set.seed(0)
figA5_B <- filter(study2_alt, condition == "COVID-19") %>% 
  ggplot(aes(x = vignette_similarity, y = allocation_normalized, fill = condition)) +
  geom_jitter(col = palette_to_use[3], alpha = 0.3, width = 0.25, height = 0.25) +
  geom_smooth(formula = y ~ x, method='lm', color = "black", size = 0.75, fill = palette_to_use[3], alpha = 1) +
  scale_x_continuous(breaks = seq(0,10,1)) +
  labs(x = "", y = "", title = "(B) COVID-19 similarity") +
  theme + 
  theme(strip.background = element_blank()) + 
  scale_y_continuous(labels = scales::percent) +
  coord_cartesian(ylim = c(0, 3.5)) + 
  guides(fill = FALSE, color = FALSE)

set.seed(0)
figA5_C<- filter(study2_alt, condition == "COVID-19") %>% 
  ggplot(aes(x = vignette_deservingness, y = allocation_normalized, fill = condition)) +
  geom_jitter(col = palette_to_use[c(10)], alpha = 0.3, width = 0.25, height = 0.25) +
  geom_smooth(formula = y ~ x, method='lm', color = "black", size = 0.75, fill = palette_to_use[10], alpha = 1) +
  scale_x_continuous(breaks = seq(0,10,1)) +
  labs(x = "Subjective deservingness", y = "Allocation\n(% of mean, $198)", title = "(C) GST deservingness") +
  theme + 
  theme(strip.background = element_blank()) + 
  scale_y_continuous(labels = scales::percent) +
  coord_cartesian(ylim = c(0, 3.5)) + 
  guides(fill = FALSE, color = FALSE)

set.seed(0)
figA5_D <- filter(study2_alt, condition == "COVID-19") %>% 
  ggplot(aes(x = vignette_similarity, y = allocation_normalized, fill = condition)) +
  geom_jitter(col = palette_to_use[10], alpha = 0.3, width = 0.25, height = 0.25) +
  geom_smooth(formula = y ~ x, method='lm', color = "black", size = 0.75, fill = palette_to_use[10], alpha = 1) +
  scale_x_continuous(breaks = seq(0,10,1)) +
  labs(x = "Subjective similarity", y = "", title = "(D) GST similarity") +
  theme + 
  theme(strip.background = element_blank()) + 
  scale_y_continuous(labels = scales::percent) +
  coord_cartesian(ylim = c(0, 3.5)) + 
  guides(fill = FALSE, color = FALSE)


(figA5_A + figA5_B) / (figA5_C + figA5_D)

ggsave("outputs/figA5.pdf", width = 7.5, height = 4)
ggsave("outputs/figA5.png", width = 7.5, height = 4)

```

### F-test alt study

```{r}

study2_alt %>% 
  cj_anova(.,
           id = ~id, by = ~condition,
           formula = allocation_normalized ~ vignette_ethnicity + vignette_health + 
             vignette_citizen + vignette_children + vignette_employment + vignette_income)

```

### Table A4: Subjective

```{r}

tabA4_models <- list()
tabA4_models[[1]] <-  miceadds::lm.cluster(data = filter(study2_alt, condition == "COVID-19"),
                  vignette_allocation ~ vignette_deservingness + vignette_similarity + id, cluster = "id")
tabA4_models[[2]] <- miceadds::lm.cluster(data = filter(study2_alt, condition == "GST"),
                  vignette_allocation ~ vignette_deservingness + vignette_similarity + id, cluster = "id") 
tabA4_models[[3]] <- miceadds::lm.cluster(data = study2_alt,
                  allocation_normalized ~ condition*vignette_deservingness + condition*vignette_similarity + id, cluster = "id") 

texreg(
  list(tabA4_models[[1]]$lm_res, tabA4_models[[2]]$lm_res, tabA4_models[[3]]$lm_res),
  omit.coef = "id",
  file = "outputs/tabA4.tex",
  single.row = TRUE,
  custom.model.names = c("COVID-19",
                         "GST rebate",
                         "Mean-normalized"),
  custom.coef.names = c("Constant",
                        "Deservingness (0-10)",
                        "Similarity (0-10)",
                        "GST",
                        "GST x Deservingness (0-10)",
                        "GST x Similarity (0-10)"),
  reorder.coef = c(1,4,2,3,5,6),
  dcolumn = TRUE,
  use.packages = FALSE,
  stars = c(0.01),
  label = "tab:a4",
  caption = "Subjective evaluations of deservingness and similarity",
  caption.above = TRUE,
  custom.note = "\\parbox{1\\linewidth}{\\vspace{2pt}%stars. Linear regression for subjective evaluations of deservingness and similarity, with individual respondent controls and clustered standard errors in parentheses. Dependent variable: allocation of cash transfer to hypothetical individuals, either under the COVID-19 or GST conditions.}"
)
```

## Appendix E: Deservingness and income

### Figure A6: Deservingness

```{r}

figA6_A <- bind_rows(
  filter(study2, condition == "COVID-19") %>% 
    cj(., id = ~id, estimate = "amce", level_order = "ascending",
         vignette_deservingness ~ vignette_gender + vignette_ethnicity + vignette_health + vignette_citizen + vignette_marital + vignette_children + vignette_employment + vignette_income) %>% 
    mutate(condition = "COVID-19"),
  filter(study2, condition == "GST") %>% 
    cj(., id = ~id, estimate = "amce", level_order = "ascending",
         vignette_deservingness ~ vignette_gender + vignette_ethnicity + vignette_health + vignette_citizen + vignette_marital + vignette_children + vignette_employment + vignette_income) %>% 
    mutate(condition = "GST")
) %>% 
  mutate(feature = as.character(feature)) %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  filter(estimate != 0) %>% 
  refactor_amce() %>% 
  ggplot(aes(x = level, y = estimate, fill = condition, col = condition)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
  geom_point(position = position_dodge(0.8)) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, position = position_dodge(0.8)) +
  theme +
  coord_flip() +
  scale_color_manual(values = palette_to_use[c(3,10)]) + 
  scale_x_discrete(position = "top", labels = NULL) + 
  facet_grid(feature~., scales = "free", switch="both") +
  labs(x = "", y = "AMCEs (0-10 deservingness scale)", title = "(A) AMCE") +
  theme(strip.background = element_blank(), 
        axis.ticks = element_blank(),
        strip.text = element_text(size = 8),
        legend.title = element_blank(), legend.position = "bottom")


figA6_B <- study2 %>% 
  amce_diffs(., id = ~id, by = ~condition, level_order = "descending",
         vignette_deservingness ~ vignette_gender + vignette_ethnicity + vignette_health + vignette_citizen + vignette_marital + vignette_children + vignette_employment + vignette_income) %>% 
  mutate(lb = estimate - 1.96*std.error,
         ub = estimate + 1.96*std.error) %>% 
  refactor_amce() %>% 
  ggplot(aes(x = level, y = estimate, col = "Difference")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgrey") +
  geom_point(col = "black") +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0, col = "black") +
  theme +
  coord_flip() +
  scale_x_discrete(position = "top") + 
  scale_color_manual(values = "black") +
  facet_grid(feature~., scales = "free", switch="both") +
  labs(x = "", y = "AMCE difference (0-10\ndeservingness scale,\nGST minus COVID-19)",
       title = "(B) AMCE differences") +
  theme(strip.background = element_blank(), strip.text = element_blank(), axis.ticks = element_blank()) + theme(legend.position = "bottom", legend.title = element_blank())

figA6_A + figA6_B

ggsave("outputs/figA6.pdf", width = 7.5, height = 8.5)
ggsave("outputs/figA6.png", width = 7.5, height = 8.5)

```

### Table A7: Deservingness, income, GST condition

```{r}

tabA7_models <- list()

tabA7_models[[1]] <- lm(data = study2,
                  vignette_deservingness ~ condition)

tabA7_models[[2]] <- miceadds::lm.cluster(data = study2,
                  vignette_deservingness ~ condition*vignette_income + id, 
                  cluster = "id")

tabA7_models[[3]] <- miceadds::lm.cluster(data = study2,
                  allocation_normalized ~ condition*vignette_deservingness*vignette_income + id, 
                  cluster = "id")

texreg(
  list(
    tabA7_models[[1]],
    tabA7_models[[2]]$lm_res,
    tabA7_models[[3]]$lm_res
  ), 
  omit.coef = "id", 
  single.row = TRUE,
  file = "outputs/tabA7.tex",
  custom.model.names = c(
    "1: Deservingness",
    "2: Income",
    "3: Allocation"
  ),
  custom.coef.names = c(
    "Constant",
    "GST",
    "\\$90,000",
    "\\$60,000",
    "\\$30,000",
    "GST x \\$90,000",
    "GST x \\$60,000",
    "GST x \\$30,000",
    "Deservingness",
    "GST x Deservingness",
    "Deservingness x \\$90,000",
    "Deservingness x \\$60,000",
    "Deservingness x \\$30,000",
    "GST x Deservingness x \\$90,000",
    "GST x Deservingness x \\$60,000",
    "GST x Deservingness x \\$30,000"
  ),
  reorder.coef = c(1,5,4,3,
                   2,8,7,6,
                   9,10,13,12,11,
                   16,15,14),
  groups = list(
    "Recipient income (ref \\$120,000)" = 2:4,
    "GST condition (ref \\$120,000)" = 5:8,
    "Deservingness (ref \\$120,000)" = 9:13,
    "Triple interaction (ref \\$120,000)" = 14:16
  ),
  dcolumn = TRUE,
  use.packages = FALSE,
  stars = c(0.01),
  label = "tab:a7",
  caption = "Non-registered tests on the relationship between deservingness, income, and the GST condition",
  caption.above = TRUE,
  scalebox = 1,
  float.pos	= "!h",
  custom.note = "\\parbox{0.9\\linewidth}{\\vspace{2pt}%stars. Linear regressions, with individual respondent controls and clustered standard errors in parentheses (for Models 2 and 3). Dependent variable for Models 1 and 2 is subjective evaluation of deservingness. Dependent variable for Model 3 is mean-normalized allocation to hypothetical individuals, either under the COVID-19 or GST conditions.}"
)

```

```{r}

```


